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= 摘要 :， 系 外 行星 的 长 期 动力 学 演化 需要 可 靠 的 数值 计算 方法 。 辛 算法 具有 保 能 量 、 保 结构 的 特 
点 ， 是 研究 哈密 顿 系统 长 期 演化 的 最 佳 积分 工具 。 双 自转 系 外 行星 后 牛顿 哈密 顿 系统 中 坐标 和 动 
量 不 可 分 离 ， 显 式 辛 算法 不 能 直接 应 用 。 利 用 相 空 间 扩充 构造 显 式 类 辛 算法 不 会 引入 人 工 耗 散 ， 
A 4 有 保 能 量 的 优势 。 主 要 探讨 四 阶 中 点 置换 相 空间 扩充 显 式 类 辛 算法 在 双 自转 系 外 行星 后 牛顿 
© 轨道 中 的 数值 性 能 。 结 果 表明 ， 中 点 置换 相 空 间 扩充 显 式 类 辛 算法 ， 在 顺 行 非 近 共 面 轨道 的 算法 
精度 与 RKF8(9) 算法 精度 差 一 个 量 级 ;在 逆行 轨道 和 近 圆 轨道 的 算法 精度 与 RKF8(9) 算法 精 
度 相 当 ; 在 偏心 率 小 于 0.9 的 轨道 表现 出 高 于 同 阶 对 比 算 法 的 算法 精度 。 除 此 之 外 ， 该 算法 具有 
良好 的 稳定 性 ， 计 算 效率 约 是 RKF8(9) 的 3 fi. 
Al x Mil: 系 外 行星 ， 隐 式 中 点 法 ， 相 空间 扩充 类 辛 算法 
© 中 图 分 类 号 :P132 文献 标识 码 : A 


1 引 言 


= 系 外 行星 的 研究 在 1995 年 迎 来 突破 ，Mayor 和 Queloz" 公布 了 有 史 以 来 发 现 的 第 一 颗 
p 环绕 类 太阳 恒星 运行 的 系 外 行星 (51 Peg b)， 这 是 人 类 观测 到 的 第 一 颗 热 木星 。 热 木星 是 一 
类 质量 在 0.3My < m(m x sini) < 13M; (M; 表示 木星 质量 ，; 表示 轨道 面 法 向 与 观测 者 视 
线 方向 的 夹 角 )， 轨 道 周期 p < 10 d 的 气 巨 星 。 和 截止 2021 年 12 月 8 日 ， 已 经 认证 的 系 外 行 
星 达到 4576 NEL KMRL A 492 颗 。 由 系 外 行星 的 观测 数据 可 知 ， 系 外 行星 与 太阳 系 
内 行星 在 轨道 特征 上 明显 不 同 。 热 木星 、 热 地 球 和 高 偏心 率 系 外 行星 的 发 现 都 对 传统 的 行星 
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成 理论 提出 了 挑战 。 通 过 Rossiter-McLaughlin" ”效应 发 现 ， 一 些 热 木星 的 轨道 面 法 向 与 
主星 的 自转 轴 之 间 的 夹 角 很 大 ， 甚 至 存在 一 些 逆行 的 热 木星 一。 
目前 研究 热 木 星 高 倾角 形成 的 机 制 都 是 在 牛顿 力学 的 框架 内 进行 的 ， 但 是 在 相 邻 行星 
或 伴星 的 倾斜 轨道 下 ， 轨 道 - 自 旋 耦 合 后 牛顿 项 会 导致 主星 的 自 旋 进 动 。 高 阶 后 牛顿 近似 主 
要 分 为 后 牛顿 拉 格 朗 日 近似 方程 、 自 洽 方程 及 ADM 坐标 系 下 的 后 牛顿 哈密 顿 量 ””。 在 
同 阶 情况 下 ，ADM 坐标 系 下 的 后 牛顿 哈密 顿 量 和 后 牛顿 拉 格 朗 日 形式 的 物理 等 价 性 已 经 得 
到 验证 。 但 是 ， 这 两 种 后 牛顿 形式 在 轨道 动力 学 行为 上 并 不 是 完全 等 价 的 ™”™。 当 
两 种 形式 中 的 一 种 通过 勒 让 德 变换 转化 为 男 一 种 形式 时 ， 高 阶 后 牛顿 项 被 截断 。 对 于 像 太阳 
系 这 样 的 弱 引 力 场 ， 这 种 差异 可 以 忽略 不 计 ， 但 对 于 致密 天 体 这 样 的 强 引 力 场 ， 可 能 会 对 两 
种 形式 的 动力 学 产生 重要 影响 。 有 两 种 运动 方程 可 以 求解 保守 的 拉 格 朗 日 动力 学 系统 ， 分 别 
为 截断 近似 的 运动 方程 和 无 截断 自治 运动 方程 ， 这 两 者 之 间 也 存在 差异 。 这 种 差异 可 
能 使 得 两 种 后 牛顿 近似 具有 不 同 的 可 积 性 和 不 可 积 性 ， 或 有 序 性 和 混沌 性 。Wu 和 Xie" 提 
出 的 正则 共 斩 自 旋 变 量 ， 在 判定 自 旋 致 密 双 星 后 牛顿 哈密 顿 系统 的 可 积 性 或 不 可 积 性 时 起 
着 非常 重要 的 作用 ， 由 此 论证 了 后 牛顿 单 体 自 旋 致 密 双 星系 统 是 可 积 的 。 
无 论 是 在 牛顿 引力 理论 还 是 在 广义 相对 论 情形 下 ， 系 外 行星 大 偏心 率 或 大 倾角 轨道 ( 行 
星 轨 道 面 法 向 与 主星 自转 轴 之 间 的 夹 角 ) 的 长 期 演化 过 程 都 需要 可 靠 的 数值 计算 方法 。 早 期 
求解 微分 方程 的 泰勒 级 数 法 ”和 龙 库 法 在 长 期 轨道 积分 时 并 不 可 靠 。Ruth” 和 冯 康 ”分 
别 在 1983 年 和 1984 年 提出 辛 算法 ， 该 算法 基于 哈密 顿 力学 的 基本 原理 ， 使 离散 化 的 差分 
方程 保持 哈密 顿 系统 的 原始 辛 结构 ， 具 有 长 期 稳定 性 。 辛 算法 的 提出 使 数值 计算 方法 的 发 
展 进入 了 全 新 阶段 。 在 太阳 系 内 ，Wisdom-Holman 辛 算法 “最 为 常用 ， 该 算法 是 将 哈密 顿 
系统 分 解 为 主要 部 分 和 次 要 部 分 进行 计算 。Laskar 和 Robutel 发 展 了 高 阶 辛 算 法 ， 并 用 上 其 
BRB AA. Ruth 提出 的 显 式 辛 算法 只 能 用 于 求解 坐标 、 动 量 可 分 离 的 哈密 顿 系统 ， 
计算 效率 和 稳定 性 较 高 。 对 于 变量 不 可 分 离 的 哈密 顿 系统 ， 显 式 辛 算法 不 能 直接 应 用 ， 可 
FA Rash SEA ERA EET, USMS AAEM ERE. 
ET 提出 的 以 隐 式 中 点 法 为 基础 的 隐 式 辛 算法 ， 适 用 于 求解 任何 哈密 顿 系统 ， 但 是 该 算法 
需要 进行 多 次 迷 代 才能 达到 较 高 的 精度 ， 计 算 效 率 明显 低 于 显 式 辛 算法 ， 一 般 不 作为 第 一 选 
择 。 相 空间 扩充 显 式 类 辛 算法 最 早 由 Pihajoki"" 提出， 通过 使 用 动量 置换 使 得 显 式 辛 算法 
在 含有 不 可 分 变量 的 哈密 顿 系统 中 的 应 用 成 为 可 能 ， 但 是 破坏 了 哈密 顿 系统 的 辛 结构 ， 不 
属于 严格 意义 上 的 辛 算法 ， 故 称 之 为 类 辛 算法 。Tao” 在 此 基础 上 提出 将 哈密 顿 系统 分 为 三 
个 部 分 ， 保 持 了 算法 的 辛 结构 。 后 续 Liu AT" Pihajoki 的 基础 上 找到 了 更 好 的 置 
换 方式 ， 结 合 Yoshida™ 高 阶 算法 的 构建 ， 提 出 了 高 阶 坐 标 动量 置换 相 空间 扩充 显 式 类 辛 算 
法 。 该 算法 通过 两 次 置换 (坐标 置换 和 动量 置换 ) 之 后 可 以 达到 高 精度 和 稳定 性 ， 但 在 后 牛 
顿 致密 双星 系统 中 的 一 些 混沌 轨道 无 法 保持 系统 的 能 量 稳定 “™。 为 了 解决 这 个 问题 ，Luo 
等 人 ™ "提出 中 点 置换 相 空 间 扩充 显 式 类 辛 算法 ， 该 算法 只 需要 进行 一 次 中 点 置换 就 能 i 
到 较 高 的 精度 和 稳定 性 。Wu 5E AP" 在 扩大 相 空 间 的 基础 上 提出 优化 算法 ， 与 未 优化 算法 相 
比 ， 其 精度 得 到 提升 。Li 和 Wo Tft Mikkola 和 Palmer" 的 基础 上 再 结合 各 种 置换 方法 ， 
提出 相 空 间 扩充 对 数 哈密 顿 显 式 对 称 算法 ， 该 算法 在 解决 高 偏心 率 轨道 和 混沌 轨道 方面 表 
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现 E lane Pan 等 人 ”构建 的 自 洽 后 牛顿 拉 格 朗 日 方程 的 扩大 相 空间 半 


在 能 量 保持 上 也 展现 出 好 的 效果 。 近 期 ， 我 国学 者 把 黑洞 时 空 分 解 为 有 多 个 显 


式 分 析 解 的 部 分 }， 然 后 再 组 合成 显 式 辛 算法 ， 解 决 了 相对 论 黑洞 时 空 构建 显 式 辛 算法 的 国际 


难题 。 


本 文 主要 研究 中 点 置换 相 空 间 扩 充 显 式 类 辛 算法 在 双 自 转 后 牛顿 系 外 行星 轨道 中 的 应 


用 。 一 阶 后 牛顿 双 自 转 二 体 问题 哈密 顿 系统 中 ， 变 量 不 可 分 离 ， 并 且 自 旋 变 量 是 以 非 正 则 


的 形式 表示 ， 


和 Xie” 的 方法 ， 通 过 引入 类 柱 坐 标 系 实现 自 旋 变 量 正则 化 ， 使 系统 的 相 空 间 具 有 完整 的 辛 
结构 。 在 一 阶 后 牛顿 双 自 转 二 体 问 题 中 采用 四 阶 隐 式 中 点 法 、 四 阶 坐标 动量 置换 相 空 间 扩充 
AK 辛 算法 、 四 阶 中 点 置换 相 空 间 扩充 显 式 类 辛 算法 和 RKF8(9) 算法 进行 数值 模拟 。 讨 


论 四 种 算法 


不 有 具备 全 局 意义 上 的 辛 结构 ， 因 此 限制 了 辛 算法 的 进一步 应 用 。 我 们 使 用 Wu 


在 双 自 转 后 牛顿 系 外 行星 轨道 (高 倾角 轨道 、 低 俩 心率 轨道 和 高 偏心 率 轨 道 ) 中 


iis 度 和 效率 。 给 出 四 种 算法 的 能 量 误差 与 行星 轨道 倾角 和 偏心 率 的 关系 图 ， 对 算法 能 


行 了 讨论 。 


量 误差 与 轨道 倾角 或 偏心 率 的 关系 进行 讨论 ， 最 后 对 不 同 算法 数值 求解 恒星 的 视 向 速度 进 


2 一 阶 后 牛顿 双 上 自转 二 体 问 题 


2.1 ”物理 模型 


一 阶 后 和 


F 顿 双 自转 二 体 问题 简化 为 质心 参考 系 后 ， 哈 密 顿 量 可 以 写 为 一 ， 


H = Hx + Hı ; (1) 


其 中 ，e 是 微 扰 理论 的 小 参数 ，e = 1/8. ¢ 表示 光速 。 


Hy id 1J p GM 
ES ' Qu r : 


2h 


Ay 是 未 受 扰 动 的 哈密 顿 量 ，Hi 可 以 展开 为 : 


7N , 


Hı = Hiew + Hso ， (3) 
1 p GM "um DAN  @M? 
ws =n ie D 2: (34 Zr Fu Bec 2,2 (4) 
一 阶 后 牛顿 轨道 哈密 顿 量 ， 
2G [/, 3m 3m 
Hso 73 [t | sme) y, «im AE (r x p) (5) 


是 1.5PN 自 旋 -轨道 耦合 项 。 在 以 上 这 些 公 式 中 ，G 是 万 有 引力 常数 ，ml Mm 表示 两 个 


天 体 的 质量 ， 


其 中 ，m 代表 次 级 天 体 的 质量 ，mz 代表 主 天 体 的 质量 。M = mi +m 是 两 


f" desee Os 4/0 BB ll 
ChinaX iva fERBTII 


4 期 Abia, Se: 显 式 类 辛 算法 在 双 自 转 系 外 行星 后 牛顿 轨道 中 的 应 用 601 


个 天 体 的 总 质量 ,1 = mima/M,v = u/M, p = p, = —p» p, Mp, THAR DADA 
天 体 1 和 天 体 2 的 动量 。” = rm 方向 由 天 体 2 指向 天 体 1, n 表示 单位 径 撩 。 自 旋 变 量 J 
在 式 (5) 中 以 非 正则 变量 的 形式 表示 ， 因此 哈密 顿 系统 是 局 部 保 辛 的 ， 不 具有 完整 的 辛 结 
构 ， 不 利于 进一步 使 用 辛 算法 进行 数值 模拟 。Wu 和 Xie” 提 出 通过 引入 类 柱 坐 标 系 来 表示 
自 旋 变量 ， 利 用 守恒 自 旋 量 将 自 旋 变量 转化 为 关于 正则 坐标 6; 和 正则 动量 €, 的 形式 。 自 旋 
变量 正则 化 使 得 原始 哈密 顿 系统 由 12 维 空间 降 为 10 维 相 空间 ， 使 整个 哈密 顿 系统 的 相 空 
间 具 有 完整 的 辛 结构 。 自 旋 变 量 用 (pi, 0i éi) 表示 为 : 

Ji = xum?Ji , (6) 


Sih, Ji 表示 单位 自 旋 矢 量 ， 即 |J;|=1， 表 示 为 : 
pi cos 0; 
J = Pi sin 6; ; (7) 
Ki&i 
fi = l(ximi) , (8) 
pi = Vil— (mE)? : (9) 
将 式 (7) — (9) RAR (6) H, ZIFRI LES STE UHR A e ae i 
| v J? — €? cos 6; | 


J; = v J? — £? sin 0; 
E? 
因此 ， 质 心 坐标 系 下 含有 正则 自 旋 变量 的 一 阶 后 牛顿 双 自转 二 体 问 题 哈密 顿 量 表示 为 : 


#1q'q'qq(i=1,2) , (10) 


H = H(r,p, Ji (91, £1), J2(62, 2)) , (11) 
" = 其 中 ， 
-一 iJ; 1J} pP GM 
d HN L 2 2 r s) 
1 p GM p?  v(zp.- Ypy + zp) B 
Hipu ZI De |o rst ym 十 (13) 
2G 3m 
Hso = “3 {| | m) CECD «(a mu) (VB- ien) VDz 一 


(: 
zpy) + : 十 iz) (V - sino.) 十 (1+ =) (V3 -& sind) | ( 2p, — 


14 6 i (1 imi) g (ep — vpe) | . (14) 
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本 文 假设 所 有 的 天 体 都 是 围绕 自转 轴 旋 转 的 球体 ， 遵 循 Barke ”和 Wex” 的 解释 ， 
自 旋 变量 J, 可 以 近似 地 表示 为 定 轴 转 动 球体 的 旋转 角 动 量 ， 即 ; 


其 中 ，w; 是 天 体 i 的 旋转 角速度 ， 五 是 天 体 i 的 惯性 矩 。 恒 星 和 行星 的 惯性 匀 
各 自 的 质量 、 半 径 有 关 。 人 恒星 的 惯性 矩 表达 式 为 ”: 


TO] 


的 表达 式 与 


Im mpm, (16) 


B 取 值 在 0.175~0.385 Ziel. BARE EXE CAI : 


IxamR? , (17) 


当 行 星 是 气 巨星 时 ，a =0.2; 属于 固态 行星 时 ，a = 0.25。 
2.2 ”算法 描述 与 构造 

X (11) 一 (14) 中 的 坐标 与 动量 不 可 分 离 ， 显 式 辛 算法 不 能 直接 应 用 。 隐 式 辛 算法 和 相 
空间 扩充 类 辛 算法 在 含有 不 可 分 变量 的 哈密 顿 系 统 中 是 不 错 的 选择 。 常 用 的 算法 主要 有 隐 
式 中 点 法 、 坐 标 动 量 置换 相 空 间 扩充 显 式 类 辛 算法 、 中 点 置换 相 空 间 扩充 显 式 类 辛 算法 。 下 
面 主要 介绍 这 三 种 算法 ， 并 结合 测试 模型 进行 相应 的 算法 构造 。 
2.231 BAP RE 


隐 式 中 点 法 在 含有 不 可 分 离 变 量 的 哈密 顿 系统 中 的 实用 性 ， 在 很 多 算 例 中 都 得 到 了 证 
明 。 我 们 可 以 使 用 牛顿 兴 代 法 和 赛 德 尔 迭 代 法 直接 进行 求解 ， 从 第 n 步 到 第 (n 十 1) 步 的 差 
分 离散 格式 写 为 : 

| OH dn+1 T dn 0,44 "WC On Pn+1 T Pn Enti m (m 
dn+1 = dn 71 h Op ( 2 , 2 , 2 , 2 
a ONE n (Suc an Past 0, Pott Pn ma ss) 
IM2: : 18 
— ,9H Qn41 T dn On41 ae On Pn+1 T Pn E aid e Es ) 
Pari = Pa Oe n a G 
© OH Qn41 + dn 0,4 T 0,, Pnii T Pn Enti WO m 
Ss Fl en h 00 ( 2 , 2 , 2 D 2 
Yoshida ”提出 了 高 阶 算法 的 构建 ， 主 要 是 在 二 阶 显 式 辛 算法 的 基础 上 ， 利 用 三 种 2n 


阶 算法 在 合适 的 步 长 系数 组 合 下 对 称 构造 一 个 (2n + 2) 阶 算法 : 


Son42(h) = Son(Arh)San(Azh)San(Arh) ， (19) 


RBA, 和 AQ A: 
_91/(2n+1) 1 


à= 9 — 91/@n+1)’ 27 9 _ 91/@n+1) 


(20) 


a 
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通过 使 用 Yoshida ”构建 高 阶 算法 的 方式 ， 四 阶 隐 式 中 点 法 IMA 的 表达 式 可 以 写 为 : 


IMA(h) = IM204h)IEM205h)IM2(AR) ， (21) 
系数 和 和 A 为 : 
—91/3 1 
^-7g-mso "cy Q2) 


+ 


2.2.2 ” 相 空 间 扩 充 类 辛 算 法 

Pihajoki 在 2015 年 提出 了 相 空 间 扩充 类 辛 算 法 。 我 们 通过 对 原始 哈密 顿 量 进行 相 空 
间 扩 充 ， 复 制 原始 变量 得 到 (4. 0,, 05, D, &y, €»), 双 倍 化 相 空 间 变量 使 得 (q, p, 01, 05, £1, E2) 
— (q, d, p, D, 91, 01, 92, 05, &, &, €», 乌 )。 最 终 得 到 两 个 完全 相等 的 哈密 顿 量 Hh (q, 0, 
05, D, &, &), He(q, bi, 92, p, &, &) "si an d H 完全 相同 ， 它 们 组 成 一 个 新 的 哈 
密 顿 量 : 


H(q,q, 0i, 0i, P, D, &,&) = H, (q, 0, D,E) + H3(q, 0i, p, &;), (i - 1,2) H (23) 
哈密 顿 量 H. 和 Hy 是 坐标 、 动 量 可 分 离 的 哈密 顿 系统 ， 相 互 独立 且 可 积 ， 可 以 分 别 进 
行 解析 求解 。 
A, 相应 的 正则 方程 写 为 : 
dg OF dh OF dôs OH 
dt Op’ dt 8€, dt 0€; 


(24) 
dp OH, dé, = OH, dé» OH, 
dt dq’ dt 06,’ dt 06, 
Ay 相应 的 正则 方程 写 为 
dq u OH» dd; u OH» dé u OH» 
dt Op’ dt Ot' dt OE, ' 
(25) 


dp | OH; d& . OH; dé OH». 


dt dq’ dt 00,' dt 00, 
二 阶 辛 算法 可 应 用 于 哈密 顿 量 H, h 是 步 长 。 二 阶 辛 算法 可 写 为 : 
h h 
s) =m (2) mom (7) | (26) 


在 初始 时 刻 ， 原 始 坐 标 、 动 量 与 扩展 后 的 坐标 、 动 量 完全 相等 ， 但 随 考 时 间 的 累积 让， 在 
积分 过 程 中 由 于 原 变量 与 扩展 变量 的 耦合 效应 ， 两 者 会 逐渐 变 得 不 相等 。Pihajoki"” 将 二 
阶 辛 算 法 进行 了 改正 ， 通 过 引入 置换 变量 M, 解决 了 这 个 问题 ， 即 : 


So(h) = M3H» (3) Hi (5) MH, (3) Hy (5) ; (27) 


Rt, Mi 和 My 代表 动量 置换 p o p. 


ve 


Chin-c 
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2.2.3 连续 坐标 动量 置换 相 空间 扩充 类 辛 算 法 
Liu & A^ " E Pihajoki"" 的 基础 上 找到 了 更 好 的 置换 方式 ， 提 出 了 四 阶 连续 坐标 动量 
置换 相 空间 扩充 显 式 类 辛 算法 ， 表 示 为 : 


S4 = S( M h)So (oh) So (Ash) M; x So (Ash) SS (Agh)SS(A,h) Ma, (28) 


HB, Mi 表示 动量 置换 p o p,& ei, 66, Mo 表示 坐标 置换 g oq. 0, o 6, 
6; ++ fos AME AL = A2 = 1/[2(2- 29), As = 5 — BAe 

连续 坐标 动量 置换 相 空 间 扩 充 显 式 类 辛 算法 ， 需 要 进行 两 次 三 重 积 分 才能 达到 高 精度 
效果 。 在 后 牛顿 哈密 顿 问题 和 旋转 哈密 顿 问题 等 哈密 顿 系统 不 可 分 问题 中 ， 该 算法 得 到 了 很 
好 的 应 用 ， 但 在 旋转 致密 双星 后 牛顿 哈密 顿 系 统 的 一 些 混 沌 轨道 中 无 法 保持 系统 的 能 量 稳 
gen 
2.24 中 点 置换 相 空间 扩充 类 辛 算法 

Luo 5: AP* " 在 Liu AT T 的 基础 上 进行 了 调整 ， 提 出 了 中 点 置换 相 空间 扩充 显 式 
类 辛 算法 。 该 算法 不 需要 进行 两 次 置换 ， 只 需要 进行 一 次 中 点 置换 就 可 以 得 到 很 好 的 精度 和 
稳定 性 。 在 积分 过 程 中 ， 中 点 置换 使 得 原始 坐标 、 动 量 与 扩展 坐标 、 动 量 始终 相等 。 四 阶 中 
点 置换 相 空 间 扩充 类 辛 算法 AA 表示 为 : 


A4 = M ® S; (sh) SS (Ah) SS (Ah), (29) 


Hp, ABA, = à = 1/(2- 23), 4,21— 224, M 代表 中 点 置换 ， 即 : 


„4+4 d „4+4 
2 "7 9- * 
46, ~ 646, 
04 a i 0; 一 : fi. 
2. 2. 
p 2 quem. 
+3 " 0) 
PtP ~ PHP 
Pp 一 2 à P 一 5 » 
& +é fog & +é 
& ? 2 & ? E 
£5 + £9 E2 + £9 
E2 > 5 ， 6 一 5 . 


3 算法 在 一 阶 后 牛顿 双 自 转 二 体 问题 中 的 应 用 


本 章 讨论 了 各 种 算法 在 一 阶 后 牛顿 双 自 转 二 体 问 题 中 的 数值 性 能 。 主 要 选取 三 类 系 外 
行星 轨道 (高 倾角 轨道 、 低 偏心 率 轨道 、 高 偏心 率 轨道 ) 进行 数值 模拟 。 我 们 分 别 采 用 四 阶 
隐 式 中 点 法 (IM4)， 四 阶 坐标 动量 置换 相 空 间 扩充 显 式 类 辛 算法 (S4)， 四 阶 中 点 置换 相 空间 
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扩充 显 式 类 辛 算法 (A4) 和 RKF8(9) 算法 在 三 类 轨道 中 进行 数值 模拟 。 将 RKF8(9) 算法 作 

为 精确 解 与 其 他 算法 进行 对 比分 析 ; 扫描 得 到 四 种 算法 的 能 量 误差 与 行星 轨道 倾角 和 偏心 

率 的 关系 图 ， 并 进行 讨论 ;最 后 对 比 了 不 同 算法 数值 求解 恒星 视 向 速度 的 精度 。 

3.1 ”初始 坐标 动量 的 获得 
初始 时 刻 的 坐标 与 动量 可 以 通过 给 定 初始 轨道 根 数 求 出 。 在 二 体 问题 中 有 5 个 不 变 的 

轨道 根 数 ， 轨 道 半 长 径 元 偏心 率 e 轨道 倾角 i、 升 交点 经 度 Q 和 近 点 角 距 ww。 还 有 轨道 

根 数 平 近 点 角 M， 会 随 着 时 间 线 性 变化 ™。 平 近 点 角 M 与 偏 近 点 角 马 满足 开 普 勒 方程 : 


E-esinE=M . (31) 
二 体 问题 的 开 普 勒 解 表示 为 : 


q = a(cos E — e)P + @vV/1 — e? sin EQ 
~2 2 


CUM sin EP + A T — e?cos EQ 
T 


3 


q= 
r 


其 中 , n= n/d 表示 平均 角 速 率 ，r 是 天 体 2 到 天 体 ZE 


om 


离 。 r 表示 为 : 


六 一 QL 一 ecos 五 ) , (32) 


拉 普 拉 斯 矢量 P 表示 为 : 


cos f2 cosw — sin f2sin w cos i 
P = | sin f2cosc + cos f2sinw cosi ; (33) 


sin w sin 


拉 普 拉 斯 矢量 Q 表示 为 : 


— cos 2sinw — sin f2 cos w cos i 
Q= | —sin Qsinw + cos f2 cos w cos i ; (34) 


cosw Sin 2 


3.2 ”数值 检验 
3.2.1 高 倾角 轨道 

轨道 1 是 高 倾角 轨道 ， 相 关 参 数 见 表 1。 其 中 M; 表示 木星 质量 ，Me 表示 太阳 质量 
Ro 表示 地 球 半 径 ，Ro 表示 太阳 半径 。 以 XO-3 b" Jgffl, XO-3 b 是 一 颗 高 倾角 热 木 星 ， 
行星 轨道 角 动 量 方向 与 主星 自转 轴 的 夹 角 投影 达到 37.3° + 3.7°. 

我 们 给 定 初始 轨道 根 数 为 区 = 0.04539 AU, e= 0.05, i= 37°, 2 = 0°, w = 0°, M = 0° 
01 = 90°, 0, = 95*， 和 初始 轨道 运动 周期 = 0.0081a， 计 算 步 长 h 设置 为 初始 轨道 周期 的 
1/600， 积 分 时 间 t = 10"T。 图 1 通过 绝对 能 量 误差 AH 展示 了 IMA, A4, S4 和 RKF8(9) 
的 能 量 误差 对 比 结果 。 信 万 = H(t) - H(0), H(t) 与 H(0) 分 别 表示 哈密 顿 量 在 时 间 上 和 0 


202306.00398v1 
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a 
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时 的 值 。 图 1 结果 显示 ， 算 法 精度 由 低 到 高 依次 为 IM4, S4 A4, RKF8(9)。 其 中 ， 扩 充 相 
空间 类 辛 算法 AA 和 SA 具有 长 期 的 稳定 性 ，A4 精度 比 SA 高 约 1 个 量 级 。S4 算法 与 IM4 
算法 精度 在 积分 前 期 几乎 一 致 ， 但 是 IM4 算法 在 计算 到 时 间 t = 106T 之 后 ， 由 于 迭代 次 数 
过 多 、 舍 入 误差 增 大 ， 因 此 能 量 误差 曲线 开始 上 移 。RKF8(9) 精度 最 高 ， 但 是 能 量 误 差 不 
稳定 并 呈 长 期 增长 趋势 ， 能 量 误差 逐渐 允 近 AA 算法 ， 积 分 时 间 到 1077 时 ， 能 量 误差 大 约 
比 A4 算法 高 1 个 量 级 。 


i ut 
表 1 ”轨道 1 的 相关 参数 从 | 
参数 WES ”数值 ”单位 ki Hoo 
mi TEME 11.7 My lli | 
ma ”主星 质量 141 Mọ Em | | | 
Rı 行星 半径 13.64 Re E | 
Ro ”主星 半 经 1.68 Ro ; 
à “轨道 半 长 径 0.04539 AU 
e 偏心 率 0.05 
i 轨道 倾角 37 (°) 
Q ” 升 交 点 经 度 0 (5 
A 行星 自 旋 倾角 i (°) 72 0 2 4 6 
^ 主星 自 旋 倾 角 1 — (9) lg(t/T) 


1 ”轨道 1 的 算法 能 量 误差 对 比 


表 2 ”轨道 2 的 相关 参数 3.2.2 ” 低 偏 心率 轨道 


参数 物理 含 数值 单位 轨道 2 是 低 偏 心率 轨道 ， 相 关 参 数 见 表 2. 

mi 行星 质量 046 M, Y M; 表示 木星 质量 ，Mo 表示 太阳 质量 ，Re 表 
ma 星 质量 1.05 Mo 示 太 阳 半 径 。 以 51 Peg p" Jl, 51 Peg b 是 一 颗 
Ri 行星 半径 0.00042 — Ro 热 木星 ， 偏 心率 仅 为 0.004 2. 

us ep. ope. MO RAIA s ML HG YUM @ = 0.05235 AU, 
ü 轨道 半 长 径 ^ 0.05235 AU f : Á N , 
nes HA e = 00042, i = 0, Q = 0, w = 0, M = 0°, 
i 轨道 倾角 0 Qj 426,0, 初始 轨道 运动 周期 了 = 0.0169 a, 
Q — 升 交点 经 度 0 (°) 计算 步 长 h 设置 为 初始 轨道 周期 的 1/400, $44) 
M ”行星 自 旋 倾 角 1 (°) 时 间 上 = 10°T. Al 2 通过 绝对 能 量 误 差 AH 展示 
A. ”主星 自 旋 倾 角 1 (°) T IM4, A4, S4 和 RKF8(9) 的 能 量 误差 对 比 结果 。 

Prot 主星 自 旋 周期 30 d 


AH = H(t) — H(0), H(t) &VIH (0) 分 别 表 示 哈 密 顿 

量 在 时 间 上 和 0 时 的 值 。 从 图 2 中 可 以 看 出 ， 算 法 
精度 由 低 到 高 依次 为 S4, IM4, A4, RKF8(9)。 扩 充 相 空 间 A4 和 SA 两 种 算法 能 量 误差 稳定 
且 没 有 发 生 误差 偏 移 ， 其 中 A4 算法 精度 比 SA 算法 高 约 两 个 量 级 。IM4 算法 能 量 误差 在 计 
算 到 时 间 上 = 10*T 时 曲线 开始 上 移 ， 误 差 变 大 。RKF8(9) 精度 最 高 ， 但 是 能 量 误差 稳定 性 


A. 
DH 
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较 差 ， 呈 直线 上 升 的 趋势 并 且 产 生 较 大 的 误差 偏 移 ， 在 积分 时 间 内 最 终 能 量 误 差 与 A4 算法 
一 样 。 
3.2.3 高 偏心 率 轨道 

轨道 3 是 高 偏心 率 轨道 ， 相 关 参 数 见 表 3。 其 中 My 表示 木星 质量 ，M。 表示 太阳 质量 ， 
Re 表示 地 球 半径 ，Ro 表示 太阳 半径 。 以 Kepler-1656 b^ 为 例 ，Kepler-1656 b 是 一 颗 高 
偏心 率 亚 土星 ， 偏 心率 达到 0.84 士 0.1。 


-8 
表 3 轨道 3 的 相关 参数 

-10 yi LT 参数 物理 含义 ”数值 单位 
M eae mi ”行星 质量 015 M; 
号 ma ”主星 质量 1.014 Mo 
m 127 Ri 行星 半径 5.02 Re 
Rə ”主星 半 经 105 Ro 
ENG à ”轨道 半 长 径 0197 AU 

A4 e 偏心 率 0.5 
RETEG) i 轨道 倾角 0  (?) 
-46 Q 升 交点 经 度 0 (°) 
-2 0 2 4 6 A 行星 自 旋 倾 角 1 (9?) 
lg(t/T) A. 主星 自 旋 倾角 1 (°) 


2 ”轨道 2 的 算法 能 量 误差 对 比 


我 们 给 定 初始 轨道 根 数 为 立 = 0.197 AU, e = 0.5, i = 0°, 2 = 0°, w = 0°, M = 0°, 
01 = 0 = 0"， 初 始 轨道 运动 周期 为 了 = 0.0868a， 计 算 步 长 h 设置 为 轨道 周期 的 1/2700， 
积分 时 间 t = 1057。 图 3 通过 绝对 能 量 误差 AH 展示 了 IM4, A4, S4 和 RKFS(9) 的 能 量 
RANT aR. AH = H(t) - H(0), H(t) M H(0) 分 别 表示 哈密 顿 量 在 时 间 上 和 0 时 的 值 。 
从 图 3 中 可 以 看 出 ， 与 RKF8(9) 精度 相 比 ，IM4, A4 和 SA 的 精度 较 低 ， 三 种 算法 的 精度 
相差 不 大 。A4 算法 精度 要 略 高 于 SA 算法 和 IM4 算法 ， 并 日 表现 出 较 好 的 稳定 性 。S4 算法 
能 量 误差 变化 波动 较 大 ， 不 具有 长 期 稳定 性 。IM4 算法 在 积分 后 期 曲线 上 移 ， 误 差 增 加 。 

综 上 所 述 ， 高 阶 龙 格 库 塔 RKF8(9) 作为 对 比 算法 ， 在 三 类 系 外 行星 轨道 都 表现 出 最 好 
的 数值 精度 ， 但 是 该 算法 稳定 性 较 差 。S4 算法 在 三 类 轨道 中 精度 相对 较 差 ， 在 低 偏心 率 轨 
道 和 高 倾角 轨道 的 情况 下 S4 算法 表现 出 较 好 的 稳定 性 ， 但 是 在 高 偏心 率 轨道 情况 下 能 量 误 
差 产 生 较 大 的 波动 。IM4 算法 精度 较 差 ， 在 低 偏心 率 轨道 和 高 倾角 轨道 后 期 ， 因 为 迭代 次 
数 多 ， 售 入 误差 增 大 ， 产 生 了 误差 偏 移 。A4 算法 在 三 类 轨道 都 表现 出 较 好 的 精度 和 稳定 性 。 


d A 列 出 了 四 种 算法 求解 三 类 轨道 所 需 的 CPU 时 间 ， 可 以 看 出 AA 的 计算 效率 明显 高 于 其 
他 算法 ，S4 算法 因为 两 次 置换 导致 计算 时 间 差 不 多 是 A4 算法 的 2 倍 ，IM4 和 RKFS(9) 需 
要 花费 更 多 的 时 间 ，RKF8(9) 所 需 时 间 约 是 A4 的 3 倍 。 总 的 来 说 ，A4 算法 不 论 是 在 精度 
方面 还 是 稳定 性 和 计算 效率 方面 都 表 现 最 好 。 
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‘A vu | 


| 


| 
表 4 四 类 算法 求解 三 类 轨道 所 需 的 
CPU 时 间 8 
轨道 A4 S4 IM4 RKF8(9) 
1 5394 9963 22300 14764 
2 4209 8900 16050 14916 
3 3425 6227 9127 9874 


lg|AH| 


lg(t/T) 
3 ”轨道 3 的 算法 能 量 误差 对 比 


3.2.4 行星 轨道 倾角 或 偏心 率 与 算法 能 量 误差 的 关系 

为 了 更 加 清楚 地 了 解 行 星 轨道 倾角 或 偏心 率 与 算法 能 量 误差 的 关系 ， 我 们 扫描 得 到 
AA 算法 、S4 算法 、IM4 算法 和 RKF8(9) 算法 的 能 量 误差 与 行星 轨道 倾角 或 偏心 率 的 关系 
图 。 在 算法 能 量 误差 与 行星 轨道 倾角 的 关系 中 ， 我 们 给 定 初始 轨道 根 数 为 立 = 0.04539 AU, 
e = 0.05, 2 = 0°, w = 0°, M = 0°, 06 = 905, 05 =95"， 初 始 轨道 运动 周期 为 卫 = 0.008 1a, 
计算 步 长 h 设置 为 轨道 周期 的 1/600， 积 分 时 间 t = 1047。 轨 道 倾 角 :的 取 值 从 O° 开始 ， 
以 8° 的 间隔 递增 到 180°. Al 4 通过 平均 能 量 误差 E 展示 了 IM4, A4, S4 和 RKF8(9) 的 能 
量 误差 与 行星 轨道 倾角 的 关系 。 从 图 中 我 们 可 以 看 出 ，RKF8(9) 的 算法 精度 几乎 不 受 行星 
轨道 倾角 的 影响 ， 算 法 精度 最 高 。 其 余 三 种 算法 ， 在 轨道 倾角 为 0"~28。 时 ， 算 法 精度 发 生 
变化 ， 大 约 提高 两 个 量 级 。 后 期 在 顺 行 非 近 共 面 轨道 ， 随 着 行星 初始 轨道 倾角 的 增加 ， 各 
算法 的 精度 趋 于 稳定 。 其 中 中 点 置换 相 空 间 扩充 类 辛 算法 AA 在 顺 行 非 近 共 面 轨道 的 精度 与 
RKF8(9) 相差 一 个 量 级 ， 而 在 逆行 轨道 中 ，A4 算法 精度 与 RKF8(9) 算法 精度 相当 。 

在 算法 能 量 误差 与 行星 偏心 率 的 关系 中 ， 我 们 给 定 初始 轨道 根 数 为 & = 0.04539AU， 
i 二 0°, 2=0°,w=0°, M = 0°, 01 = 90°, 04 = 95*， 初 始 轨 道 运动 周期 为 T= 0.0081 a, 
WEAK h 设置 为 轨道 周期 的 1/600， 积 分 时 间 t = 1047。 偏 心率 e 的 取 值 从 0 开始 ， 以 
0.1 的 间隔 递增 到 1。 图 5 通过 平均 能 量 误差 E 展示 了 IM4, A4, S4 和 RKFS(9) 的 能 量 误 
差 与 行星 偏心 率 的 关系 。 从 图 中 我 们 可 以 看 出 ， 整 体 上 RKF8(9) 算法 的 精度 最 高 。 其 余 三 
种 算法 ， 在 偏心 率 小 于 0.9 的 情况 下 ，A4 算法 具有 最 好 的 算法 精度 。 偏 心率 为 0.9 的 时 候 ， 
AA 算法 精度 受到 影响 ， 算 法 精度 由 高 到 低 依次 为 S4, IMA, A4。 偏 心率 在 0~0.1 的 时 候 ， 
IM4 算法 精度 高 于 SA 算法 ， 但 是 后 期 随 着 偏心 率 的 越 来 越 大 ，S4 算法 和 IM4 算法 的 算法 
精度 几乎 一 致 。 在 高 偏心 率 轨道 ，S4, IM4, A4 算法 精度 明显 低 于 RKF8(9) 算法 ; 这 是 因 
为 S4, IM4, A4 算法 是 定 步 长 算法 ， 在 导数 变化 剧烈 的 时 候 ， 定 步 长 会 引入 大 量 误 差 ， 而 
RKF8(9) 是 变 步 长 算法 ， 可 以 自行 调节 步 长 减少 误差 累积 。 
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图 4 ”四 种 算法 的 能 量 误差 与 行星 轨道 倾角 的 关系 图 5 ”四 种 算法 的 能 量 误差 与 行星 偏心 率 的 关系 


3.25 各 算法 得 到 的 恒星 视 向 速度 

我 们 使 用 A4, S4, IM4 和 RKF8(9) 算法 分 别 在 三 类 系 外 行星 轨道 (高 倾角 轨道 、 低 偏心 
率 轨 道 和 高 偏心 率 轨 道 ) 中 进行 模拟 ， 得 到 各 算法 对 恒星 视 向 速度 的 影响 。 图 6 是 积分 10 
个 周期 的 对 比 图 ， 纵 坐标 表示 各 算法 得 到 的 沿 轴 x 的 恒星 视 向 速度 。 从 图 中 可 以 看 出 ， 在 
三 类 轨道 中 ， 各 算法 得 到 的 沿 轴 z 的 恒星 视 向 速度 曲线 重合 。 


v. / (AU:a) 


t/T 
c) 
图 6 ”四 种 算法 得 到 的 沿 轴 x 的 恒星 视 向 速度 


后 续 进 行 更 长 时 间 的 积分 ， 以 RKF8(9) 算法 得 到 的 沿 轴 zx 的 恒星 视 向 速度 作为 精确 
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解 ， 主 要 讨论 积分 后 (高 倾角 轨道 积分 1057， 低 偏心 率 轨道 积分 107?T， 高 偏心 率 轨道 积分 
10°T) 第 一 个 周期 A4 算法 、S4 算法 、IM4 算法 得 到 的 沿 轴 x 的 恒星 视 向 速度 的 最 大 值 的 相 
对 误差 。 结 果 见 表 5。 从 表 5 中 的 数据 可 看 出 ， 三 种 算法 得 到 的 沿 轴 z 的 恒星 视 向 速度 最 大 
值 的 相对 误差 都 很 小 ， 没 有 显著 性 差异 。 


表 5 不同 算法 得 到 的 沿 轴 x 的 恒星 视 向 速度 的 最 大 值 的 相对 误差 
轨道 A4 S4 IM4 
高 倾角 轨道 ^— 1077x107? 7.504 x 1075 一 3.49 x 1078 
低 偏 心率 轨道 5.7x10-?3 4.512x107f 511x107? 
高 偏心 率 轨道 ”6.08x10-2 3.83x10-3 一 2.659 x 107? 


A 总 结 与 展望 


本 文 主要 探究 显 式 类 辛 算法 在 双 自 转 系 外 行星 后 牛顿 轨道 中 的 应 用 。 选 取 三 类 轨道 讨 
论 了 四 阶 隐 式 中 点 法 、 四 阶 坐 标 动量 置换 相 空间 扩充 类 注 算 法 、 四 阶 中 点 置换 相 空 间 扩充 类 
辛 算法 和 RKFS(9) 在 一 阶 后 牛顿 双 自 转 二 体 问题 中 的 算法 精度 和 计算 效率 ， 并 给 出 轨道 倾 
角 或 偏心 率 与 各 算法 能 量 误差 的 关系 。 数 值 模拟 结果 表明 ， 在 三 类 轨道 中 ， 四 阶 隐 式 中 点 法 
因为 需要 进行 多 次 迭代 从 而 导致 了 计算 成 本 的 增加 ， 并 且 四 阶 隐 式 中 点 法 只 有 在 低 偏心 率 
轨道 表现 出 高 于 四 阶 坐 标 动 量 置换 相 空 间 扩充 类 辛 算法 的 精度 ， 在 高 偏心 率 轨道 和 高 倾角 
轨道 的 情况 下 精度 都 相对 较 低 ， 而 且 在 积分 后 期 ， 因 为 舍 入 误差 过 大 的 原因 导致 能 量 误差 发 
生 偏 移 。 四 阶 坐 标 动 量 置换 相 空 间 扩充 类 辛 算法 在 三 类 轨道 中 表现 出 相 较 一 般 的 精度 和 稳 
定性 ， 在 高 偏心 率 轨道 能 量 误 差 产 生 较 大 的 波动 ， 并 且 因 为 需要 两 次 置换 ， 计 算 成 本 增加 。 
四 阶 中 点 置换 相 空间 扩充 类 辛 算法 整体 上 都 表现 出 较 好 的 精度 和 稳定 性 ， 计 算 效率 最 高 。 除 
此 之 外 ， 四 阶 中 点 置换 相 空间 扩充 类 辛 算法 在 逆行 轨道 表现 优异 ， 算 法 精度 和 RKF8(9) 精 
度 相 当 ， 在 偏心 率 小 于 0.9 的 轨道 表现 出 高 于 同 阶 对 比 算法 的 精度 。 

四 阶 中 点 置换 相 空间 扩充 类 辛 算 法 在 非 近 共 面 轨道 、 逆 行 轨道 和 近 圆 轨道 都 表现 出 最 
好 的 适用 性 ， 后 续 我 们 将 采用 四 阶 中 点 置换 相 空 间 扩充 类 驻 算 法 对 热 木星 高 倾角 的 形成 进 
行 探究 。 系 外 行星 系统 的 轨道 构 型 多 种 多 样 ， 未 来 还 可 以 将 算法 应 用 于 解决 更 多 的 系 外 行 
星 轨 道 的 长 期 积分 问题 ， 使 我 们 对 系 外 行星 系统 的 动力 学 形成 和 演化 有 更 加 清晰 的 了 解 与 
认识 
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Application of an Explicit Symplectic-like Integrator in 
Post-Newtonian Orbits of Double-rotating Exoplanet System 


ZHENG Jing-jing?, WANG Ying!?, LIU Fu-yao', SUN Wei!?, WANG Ya-ru!?, 
XIAO Qian-qian!?, CHEN Fen!? 


(1. Shanghai University of Engineering Science, School of Mathematics, Physics and Statistics, Shanghai 
201620, China; 2. Shanghai University of Engineering Science, Center of Application and Research of 
Computational Physics, Shanghai 201620, China) 


Abstract: Long-term dynamical evolution of exoplanets requires reliable numerical com- 
putational methods. Symplectic algorithm has the advantage of conserving energy and 
structure, which is considered as the best integrator to simulate long-term evolution of 
Hamiltonian systems. Post-Newtonian Hamiltonian of double-rotating exoplanet system 
contains inseparable forms of coordinates and momenta, where explicit symplectic integra- 
tors cannot be applied directly. The use of phase space extension to construct the explicit 


symplectic-like algorithm does not introduce artificial dissipation and has the advantage of 
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energy conservation. This paper mainly discusses the numerical performance of the fourth- 
order explicit extended phase space symplect-like method with midpoint permutations in 
post-Newtonian orbits of double-rotating exoplanet system. The results show that the accu- 
racies of the fourth-order explicit extended phase space symplect-like method with midpoint 
permutations in prograde non near coplanar orbits are only one order of magnitude lower 
than that of RKF8(9), the accuracies of the algorithm in retrograde orbits and near circular 
orbits are equivalent to these of RKF8(9), and the accuracy of the algorithm in orbits with 
eccentricity less than 0.9 are higher than these of the same order comparison algorithms. In 
addition, the algorithm has good stability and the computational efficiency is about three 
times that of RKF8(9). 


Key words: exoplanet; implicit symplectic algorithm; extended phase space symplect-like 
method 


